!------------------------------------------------------------------------------
! TJU/Department of Mechanics, Fluid Mechanics, Code START
!------------------------------------------------------------------------------
!> @Author: Jianxin Liu shookware@tju.edu.cn
!> @Date: 2024-08-13 21:38:30
!> @LastEditTime: 2024-08-13 21:38:32
!> @LastEditors: Jianxin Liu shookware@tju.edu.cn
!> @Description:
!> @FilePath: /START_NS_3D/src/splitter/module_splitter_van_leer.f90
!> @Copyright (c) 2024 by Jianxin Liu email: shookware@tju.edu.cn, All Rights Reserved.
!------------------------------------------------------------------------------
submodule(mod_splitter) submod_splitter_van_leer
   implicit none

contains

   module subroutine SplitVanLeer(this, q_in, q_l, q_r, chi_x, chi_y, chi_z, chi_face)
      class(type_splitter_van_leer), intent(inout) :: this
      class(type_vector_base), intent(in) :: q_in(:)
      class(type_vector_base), intent(inout) :: q_l(:), q_r(:)
      class(type_vector_base), intent(in) :: chi_x, chi_y, chi_z, chi_face

      real, parameter :: eps = 1.0d-6

      integer :: n, l

      class(type_self_vector), pointer :: q_in_ptr(:), q_l_ptr(:), q_r_ptr(:)
      class(type_self_vector), pointer :: chi_x_ptr, chi_y_ptr, chi_z_ptr, chi_face_ptr

      call vector_2_selfvector(q_in, q_in_ptr)
      call vector_2_selfvector(q_l, q_l_ptr)
      call vector_2_selfvector(q_r, q_r_ptr)
      call vector_2_selfvector(chi_x, chi_x_ptr)
      call vector_2_selfvector(chi_y, chi_y_ptr)
      call vector_2_selfvector(chi_z, chi_z_ptr)
      call vector_2_selfvector(chi_face, chi_face_ptr)

      call SplitVanLeer2(this, q_in_ptr, q_l_ptr, q_r_ptr, chi_x_ptr, chi_y_ptr, chi_z_ptr, chi_face_ptr)

   end subroutine SplitVanLeer

   module subroutine SplitVanLeer2(this, q_in, q_l, q_r, chi_x, chi_y, chi_z, chi_face)
      class(type_splitter_van_leer), intent(inout) :: this
      type(type_self_vector), intent(in) :: q_in(:)
      type(type_self_vector), intent(inout) :: q_l(:), q_r(:)
      type(type_self_vector), intent(in) :: chi_x, chi_y, chi_z, chi_face

      type(type_array) :: v_in(size(q_in)), v_l(size(q_in)), v_r(size(q_in))
      type(type_array) :: k_x, k_y, k_z, k_face
      real(dp) :: vv_in(size(q_in)), vv_l(size(q_in)), vv_r(size(q_in))

      integer :: l, lsta, lend, i

      do l = 1, size(q_in)
         call q_in(l)%get_array_read(v_in(l)%p)
         call q_l(l)%get_array(v_l(l)%p)
         call q_r(l)%get_array(v_r(l)%p)
      end do

      call chi_x%get_array_read(k_x%p)
      call chi_y%get_array_read(k_y%p)
      call chi_z%get_array_read(k_z%p)
      call chi_face%get_array_read(k_face%p)

      do concurrent(l=lbound(v_in(1)%p, dim=1):ubound(v_in(1)%p, dim=1))
         do concurrent(i=1:size(q_in))
            vv_in(i) = v_in(i)%p(l)
         end do
         call SplitPoint(this, vv_in, vv_l, vv_r, k_x%p(l), k_y%p(l), k_z%p(l), k_face%p(l))
         do concurrent(i=1:size(q_in))
            v_l(i)%p(l) = vv_l(i)
            v_r(i)%p(l) = vv_r(i)
         end do
      end do

      call chi_x%restore_array_read(k_x%p)
      call chi_y%restore_array_read(k_y%p)
      call chi_z%restore_array_read(k_z%p)
      call chi_face%restore_array_read(k_face%p)
      do l = 1, size(q_in)
         call q_in(l)%restore_array_read(v_in(l)%p)
         call q_l(l)%restore_array(v_l(l)%p)
         call q_r(l)%restore_array(v_r(l)%p)
      end do

   end subroutine SplitVanLeer2

   pure subroutine SplitPoint(this, v_in, v_l, v_r, k_x, k_y, k_z, k_face)
      implicit none
      class(type_splitter_van_leer), intent(in) :: this
      real(dp), intent(in) :: v_in(:)
      real(dp), intent(out) :: v_l(:), v_r(:)
      real(dp), intent(in) :: k_x, k_y, k_z, k_face

      real(dp) :: rho, u, v, w, T
      real(dp) :: ubar, c
      real(dp) :: Ma
      real(dp) :: sign, f_mass, f_energy

      real(dp), parameter :: eps = 1.0d-6

      rho = v_in(1)
      u = v_in(2); v = v_in(3); w = v_in(4)
      T = v_in(5)

      ubar = k_x*u + k_y*v + k_z*w
      c = sqrt(T)*this%eos%ref%flg_uref

      Ma = ubar/c

      if (Ma >= 1.0d0) then
         v_l(:) = v_in(:)
         v_r(:) = 0.0d0
      elseif (Ma <= -1.0d0) then
         v_l(:) = 0.0d0
         v_r(:) = v_in(:)
      elseif (Ma > -1.0d0 .and. Ma < 1.0d0) then
         #: for s, l in list(zip(["+", "-"], ["l", "r"]))
         sign = ${s}$1.0d0
         f_mass = sign*(0.25d0*rho*c)*(Ma + sign*1.0d0)**2
         f_energy = (1.0d0 - this%eos%ref%gamma)*ubar**2 + sign*2.0d0*(this%eos%ref%gamma - 1.0d0)*ubar*c + 2.0d0*c**2
         f_energy = f_energy/(this%eos%ref%gamma**2 - 1.0d0)
         f_energy = f_energy + 0.5d0*(u**2 + v**2 + w**2)
         f_energy = f_mass*f_energy
         v_${l}$ (1) = f_mass
         v_${l}$ (2) = f_mass*(k_x*(-ubar + sign*2.0d0*c)/this%eos%ref%gamma + u)
         v_${l}$ (3) = f_mass*(k_y*(-ubar + sign*2.0d0*c)/this%eos%ref%gamma + v)
         v_${l}$ (4) = f_mass*(k_z*(-ubar + sign*2.0d0*c)/this%eos%ref%gamma + w)
         v_${l}$ (5) = f_energy
         v_${l}$ (:) = v_${l}$ (:)*k_face
         #: endfor
      end if

   end subroutine SplitPoint

end submodule submod_splitter_van_leer
